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ABSTRACT 

We present multidimensional simulations of the post-explosion hydrodynamics in three different 
15 M Q supernova models with zero, \Q~ a Zq, and Zq metallicities. We follow the growth of the 
Raylcigh- Taylor instability that mixes together the stellar layers in the wake of the explosion. Models 
are initialized with spherically symmetric explosions and perturbations are seeded by the grid. Calcu- 
lations arc performed in two-dimensional axisymmetric and three-dimensional Cartesian coordinates 
using the new Eulerian hydrodynamics code, CASTRO. We find as in previous work, that Rayleigh- 
Taylor perturbations initially grow faster in 3D than in 2D. As the Rayleigh- Taylor fingers interact 
with one another, mixing proceeds to a greater degree in 3D than in 2D, reducing the local Atwood 
number and slowing the growth rate of the instability in 3D relative to 2D. By the time mixing has 
stopped, the width of the mixed region is similar in 2D and 3D simulations provided the Rayleigh- 
Taylor fingers show significant interaction. Our results imply that 2D simulations of light curves and 
nucleosynthesis in supernovae (SNe) that die as red giants may capture the features of an initially 
spherically symmetric explosion in far less computational time than required by a full 3D simulation. 
However, capturing large departures from spherical symmetry requires a significantly perturbed explo- 
sion. Large scale asymmetries cannot develop through an inverse casacde of merging Rayleigh- Taylor 
structures; they must arise from asymmetries in the initial explosion. 
Subject headings: 



1. INTRODUCTION 

Supernova 1987A has furnished astronomers with per- 
haps the best opportunity to study a core-collapse ex- 
plosion to date. One of the most exciting discoveries to 
emerge from this event is the evidence for large-scale, 
extensive mixing in the supernova ejecta. 7— ray lines 
emitted by the decay of 56 Co were detected half a year 
earlier than ex pected for a sphe rically symmetric ex- 
plosion model (|Matz et allll988t ). Modeling the light 
curve in ID requires a large amount of 56 Ni to be mixed 
outward and H and He to be mixed inward (jUtrobinl 
l2004h . The Bochum event, in which fine-structure Ha 
lines are observed 2 weeks after the explosion, is ex- 
plained by the ejection of 10~ 3 M Q of 56 Ni moving with 
a velocity in exc ess of 4000 km s -1 in to the far hemi- 
sphere of 1987A (jHanuschik et al.lll988| ). The Rayleigh- 
Taylor (RT) instability was posited as a possible explana- 
tion for the large-scale mixing implied for this explosion. 
Many groups using both SPH and grid-based codes sim- 
ulated the evolution of 87A-likc progenitor models, first 
in 2D, then in 3D. It soon became apparent that more 
than just the Rayleigh- Taylor instability is needed to ex- 
plain the high 56 Ni velocities as well as the large degree 
of outward mixing of heavy elements and inward mix- 
ing of lighter elements. Simulations that posit modest 
initial perturbations ar e unable to rep l icate these fea- 
tures of the explosions. iKifonidis et al.l ([2006) find that 
by following the explosion from early times, and using 
large, low-order perturbations in the inner layers, they 
are able to replicate the high 56 Ni velocities and large 
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amounts of mixing seen in the explosio n. The stand 
ing accretion shock instability (SASI) dBlondin et al 
2003L iBlondin fc Mezzacappal [2006: Sch eck et all 1200 
Burrows et al.1 12006L l2007t iMarek fc Jankal l2009h pro 



vides a mechanism by which the inner parts of the ejecta 
could be deformed. 

Most studies of the Rayleigh- Taylor instability in core 
collapse supernovae use 1987A-like progenitor models. 
However, supernovae in the high-redshift universe, as 
well as many modern day supernovae which die as 
red, not blue, giants, explode with presupernova struc- 
tures which differ significantly from 1987A. This may 
greatly influence the evolution of the Rayleigh- Taylor in- 
stability in their outer lay ers. Simulations of rotating 
zcro-mctallicity supernovae (Ekstrom et al. 2008; Hirschi 
120071) indicate that they may die as extended red, rather 
than compact blue, giants. Rotation mixes CNO ele- 
ments to the base of the hydrogen-burning shell, where 
they catalyze the CNO cycle, greatly boosting the burn- 
ing rate of H. This leads to the formation of a large con- 
vective envelope, effectively turning a compact blue star 
into a giant red one. Models indicate that the shell boost 
induced convection effectively mixes together the helium 
shell and the hydrogen envelope. The dense helium shell 
has been found to be important t o the post-bounce evo - 
lution of 1987A-like progenitors ((Hammer et al.l |2009[) . 
Rotation has less of an effect on early stars with very 
low, but not zero, metallicity, s uch as the Z = 10~ 4 
Zq progenitor models studied in iJoggerst et al.l (|2010f ) 
as well as the current work. These stars remain blue. 
Solar metallicity stars end their lives as red giants, but 
retain their helium shells. Red stars should evolve differ- 
ently from 1987A-type models after the supernova shock 
is launched. They are larger, so the shock takes longer 
to traverse the star, and the Rayleigh- Taylo r instabilities 
have more time to develop l|Chevalierl 1 1 9 76[) . 
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Previous sim ulations of core-collapse explosions in 
two dimensions (lArnett et alJll989t Hachisu et alJll990t 



Frvxell et all 119911: IHerant , fc Bend 119911: iMueller et al.l 
199H IHerant fe Benzj|1992l:lHacliisu et alll992l ) havepri- 



marily focused on 87A-like progenitor models. Studies of 
Raylcigh- Taylor-induced mixing in three dimensions has 
exclusively focused on attempts to replicate observations 
of 1987A, and have found significant differences between 
simulations performed in two and three d imensions. The 
most succes sful of these 2D attempts, Kifonidi s et al.l 
(|2003l 120061) found that to reproduce the high 56 Ni ve- 
locities observed in 1987A large initial asymmetries of 
low mode order were necessary. In studies of mixing 
in both red and blue s upergiants (|Joggerst et al.l 120091 : 
IHerant fc Wooslevlfl99l mixing is found to be more vig- 
orous in red stars than in blue stars. Uoggerst et al.l 
(|2010D evolve 36 models spanning three rotation rates, 
three explosion energies, three masses, and two metal- 
licities in two dimensions, finding that mixing is more 
vigorous and goes on longer in stars that die as red as 
opposed to blue giants, in higher energy explosions, and 
in less massive stars. The l iterat ure on 3D simulations 
is more sparse. iKane et al.1 (|2000| ) find that for a single- 
mode perturbation of 10% amplitude, Raylcigh- Taylor 
fingers grow 30 — 40% faster in 3D than in 2D. The 
most rece nt paper on the post -explosion hydrodynamics 
in 1987A (lHammer et al.ll2009D use a successful explosion 
model from iScheckl (1200717 " to initialize a 3D simulation, 
along with several 2D simulations derived from merid- 
ional slices through the initial explosion model. They 
find that for this model, which is highly asymmetric in 
the inner regions, the Raylcigh- Taylor instability grows 
faster and ballistically moving clumps of 56 Ni reach the 
hydrogen envelope in 3D. These clumps are effectively 
stopped by the dense He shell in the suite of 2D models 
initialized from slices through the IScheckl (I2007D model. 

This is the first paper to simulate the post-explosion 
hydrodynamics of models with a diversity of presuper- 
nova structures in three dimensions. All models have 
the same mass (15 M Q ) and explosion energy (1.2xl0 51 
ergs at infinity), but two models explode as red giants 
and one as a more compact blue star. The red giant 
with zero metallicity effectively lacks a He shell because 
of convection arising from the hydrogen shell boost. The 
Z = 10~ 4 Zq and Z= Z & models both have helium shells, 
though the former is blue and the latter red at the time 
of explosion. This diversity of models allows us to see if 
presupernova structure has an impact on the subsequent 
evolution of nearly isotropic 2D as opposed to 3D mod- 
els. These calculations set the stage for further study 
of asymmetric explosions in these diverse presupernova 
models. 

The numerical methods employed in this study are dis- 
cussed in Section § [2J and initial models and problem 
configuration are discussed in Section § [3] Results are 
presented in Section § [4j and are discussed and compared 
with previous 3D simulations in § [5] We offer some con- 
clusions in § [51 

2. NUMERICAL ALGORITHMS 

The multidimensional simulatio ns described in this pa - 
per are performed with CASTRO (lAlmgren et al.l |2010T ). 
an Eulerian adaptive mesh refinement (AMR) hydrody- 
namics code. Time integration of the hydrodynamics 



equations in CASTRO is based on an unsplit higher-order 
piecewisc parabolic method (PPM); the additional fea- 
tures are described below. All simulations described here 
used either axisymmctric coordinates in 2D or Cartesian 
coordinates in 3D. 

2.1. Equation of State 

We followed sixteen elements, from hydrogen through 
56 Ni, so that our elemental abundances could later be 
used to compute spectra. The atomic weights and 
amounts of the elements are used to calculate the mean 
molecular weight of the gas required by the equation of 
state. The abundance of 56 Ni was followed separately 
in order to calculate the energy deposited by radioactive 
decay of 56 Ni and 56 Co. 

The equation o f state in our simulations is as described 
in lJoggerst et al.l (|2010[) . It assumed complete ionization 
and included contributions from both radiation and ideal 
gas pressure: 



p = f( P ,T) i a r 4 
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where P is the pressure, a is the radiation constant, ks 
is Boltzmann's constant, T is the temperature, p is the 
density, m p is the proton mass, p, is the mean molecular 
weight, and E is the energy. 

The function /(p, T) is a measure of the contribution 
of radiation pressure to the equation of state that ac- 
counts for the fact that radiation ceases to be trapped at 
some point. It is 1 in regions where radiation pressure 
is important, i. e. where gas is optically thick, and 
in regions where radiation pressure is unimportant, i. c. 
where gas is optically thin, with a smooth transition in 
between. The function f(p,T) takes the form 
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if p > 10 9 gm cm 3 
or T < T, 
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f(T) = e 



if p < 10 9 gm cm 3 
and T > T, 



where T neg is the temperature at which contributions 
to the pressure from radiation are 100 times less than 
that contributed by ideal gas pressure: 

Zk bP 1/3 
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Radiation pressure will begin to dominate the equa- 
tion of state above T neg without some adjustment, even 
though this is unphysical, as radiation pressure is neg- 
ligible in optically thin regions. Damping the radiation 
component of the EOS with the function f(p, T) in these 
regions provides a more physical solution. 

2.2. Radioactive Decay of 56 Ni 

Energy from the radioactive decay of 56 Ni to 56 Fe was 
deposited local ly at each mesh point, in the same manner 
as described in lJoggerst et al.l (I2010D . The rate at which 
energy from the decay of 56 Ni to 56 Co was deposited in 
the grid is given by: 

dE m = X Ni X m e- XNit q(Ni). (4) 
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The decay rate of 56 Ni, X m , is 1.315 x 1CT 6 s" 1 , and the 
amount of energy released per gram of decaying 56 Ni is 
q(Ni), which we set to 2.96 x 10 16 erg g _1 . X Ni is the 
mass fraction of 56 Ni in the cell. The mass fraction of 
56 Co at a given time can be expressed in terms of the 
mass fraction of initial 56 Ni by 

X Co = - Xn \ X Nl (e- x ^ - e" A -*), (5) 

so that the energy deposition rate from 56 Co is 

dE & = - Xn \ X m (e- X ^ - e- Ac °*))A Co g(C ). (6) 

ACo — *Ni 

We assumed a decay rate Xco = 1.042 x 1CT 7 s _1 and 
an energy per gram of decaying 56 Co, q(Co), equal to 
6.4 x 10 16 erg g" 1 . 

2.3. Gravity 

Although CASTRO supports several different approaches 
to solving for self-gravity, we used the monopole approx- 
imation for gravity in the calculat ions presented here, as 
discussed in Ijoggerst et al.l (|2010f ). Because the density 
profiles in the simulations depart very little from spher- 
ical symmetry, the gravitational field constructed using 
the monopole approximation is nearly identical to that 
found by solving the full Poisson equation for the grav- 
itational potential, and using the monopole approxima- 
tion significantly reduces the computational cost of the 
calculations. 

Gravity from a point mass located at the origin was 
also included in the gravitational potential. The point 
mass represents the compact remnant left behind by 
the SN explosion. As infalling matter crosses the zero- 
gradient inner boundary near the origin, it is added to 
this point mass. 

2.4. Adaptive Mesh Refinement 

CASTRO's AMR algorithm uses a nested hierarchy of 
logically-rectangular grids with simultaneous refinement 
of the grids in both space and time. The integration al- 
gorithm proceeds recursively, advancing coarse grids in 
time, then advancing fine grids multiple steps to reach 
the same time as the coarse grid, and finally synchroniz- 
ing the data at different levels. 

The regions of refinement evolve throughout the simu- 
lation based on user-specified refinement criteria applied 
to the solution. In the simulations presented here, re- 
gridding of all grids at level £ + 1 and above occurred 
every two level £ time steps. 

Our ref inement criter ia are based on the "error esti- 
mator" of lLohnc r (1987), which is essentially the ratio of 
the second derivative to the first derivative at the point 
at which the error is evaluated. D etails of the implemen - 
tation arc discussed more fully in lAlmgren etaLl (|2010D . 
The result is a dimensionless, bounded estimator, which 
allows arbitrary variables to be subjected to the same 
criterion for error estimation. In these calculations, den- 
sity, pressure, velocity, and the abundances of 56 Ni, He, 
and O were used as refinement variables, with the el- 
emental abundances only used in a particular region if 
their abundance was greater than 10 -3 . 



3. SIMULATIONS SETUP AND EXECUTION 

The simulations presented here were carried out in two 
distinct stages. First each stellar model was evolved in 
one dimension using the code KEPLER to the point where 
the core became unstable to collapse. It was then artifi- 
cially exploded by means of a piston located at the base 
of the oxygen shell with enough energy that the explo- 
sions had 1.2 x 10 51 ergs of energy at infinity. The models 
were evolved forward in one dimension for 20 seconds for 
models zl5 and ul5, and 100 seconds for sl5, until all 
nuclear burning had ceased. 

These profiles were then mapped onto a two- 
dimensional r-z or three-dimensional Cartesian grid in 
CASTRO, and the calculation was evolved until the shock 
reached the edge of the simulation domain. The calcula- 
tion was then stopped, the domain enlarged, and the cal- 
culation restarted within the larger domain (see § [373] be- 
low). Data to fill the new regions of the enlarged domain 
were supplied from the same one-dimensional profile that 
was used to initialize the multidimensional calculation. 

This process continued until the models had evolved 
out to radii where Rayleigh- Taylor mixing ceased and 
the star was expanding essentially homologously. We 
computed three 15 M Q models, at zero, 10~ 4 Z Q , and 
solar metallicities. 

3.1. Progenitor Models 

Three different supernova models, representing three 
different metallicities, are presented in this paper. Model 
zl5 has zero initial metallicity, and represents a Popula- 
tion III (Pop III) star. Model ul5 has Z = 10~ 4 Z®, 
and model sl5 has solar metallicity All models are de- 
rived from rotating 15 M Q progenitors. The zero and low 
metallicity mod els presented in this p aper are the same 
as pre sented in Uoggerst et al.l (|2010l) and iZhang et al.l 
(|2008[ ). The solar metallicity model, model s!5. is sim - 
ilar to the model presented in lWooslev fc Hegerl (|2007|) . 
For moderate values of rotation, the amount of rotation 
was fo und to have little effe ct on its post-explosion evo- 
lution (|Joggerst et al.l [20101 ). There are, however, pro- 
found differences between rotating and no n-rotating zero 
metallicity supernovae. Several studie s (Ekstro m et al.l 
l2008tlHirschill2007l;lJoggerst et al.ll2010[ )) have found that 
a small amount of rotation in zero metallicity massive 
stars dredges up enough carbon, nitrogen, and oxygen 
from the helium burning shell to catalyze a CNO cy- 
cle at the base of the hydrogen envelope. This leads to 
rapid burning and a shell boost that triggers convection, 
mixing most of the helium and hydrogen shells together, 
and "puffing up" the outer envelope of the star. This is 
apparent in Figure [TJ Models with a small amount of 
metals (i. e. ul5 in this paper) do not show this shell 
boost behavior, and die as blue supergiants. Rotation 
also has little effect on the prcsupcrnova structure of the 
solar metallicity model sl5, which dies as a red giant. 

3.2. Initialization of Multidimensional Data 

In mapping the radial data from KEPLER onto the mul- 
tidimensional grids in CASTRO, special care was taken to 
properly resolve the key elements of the simulations: the 
shock, the elemental shells, and the 56 Ni at the center 
of the explosion. In particular, both the 56 Ni and the 
O shell were resolved with a minimum of 16 cells at 
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Fig. 1. — Shell structure for H, He, O, Si, and Ni for our initial models. Other elements were included in our simulations but are left out 
of the figure for clarity. Note that while sl5 and ul5 (top and middle figures) show an He shell, this shell is not present in model zl5 because 
of convection arising from a hydrogen shell boost. The models are 15 solar masses, but the outer layers are essentially a continuation of 
the levels at the right hand side of the plot, and have been omitted to better show the detail in the center. 

of the previous coarsest resolution, so that after each en- 
largement the base grid continued to be 128". The sim- 
ulation then had 3 levels of refinement; in most cases 
(except for ul5) the highest level of refinement was re- 
moved at this point, reducing the finest resolution by a 
factor of 2 as well. This allowed us to tailor the size of 
the domain to the size of the region of interest, without 
wasting computational resources in regions where noth- 
ing of interest was happening. Each time the domain was 
enlarged data from the new regions of the domain were 
interpolated from the original model. 

To ensure that the strategy described above did not 
affect the results of the simulations relative to analogous 
calculations performed on domains which were larger 
throughout, we ran two-dimensional simulations in which 
the final domain size was imposed from the start. The 
results were essentially the same, although more complex 
mixing structures are visible in the case with the initially 
larger domain. The width of the mixed region remains 
the same, i.e. heavy elements are mixed out the the 
same point in radius and mass, while lighter elements 
are mixed inward to the same point in the both cases. 
The velocities obtained by the elements are the same. 
This is sufficient evidence that the strategy of enlarging 
the domain throughout the calculation did not introduce 
spurious artifacts. 



the highest level of refinement. This mapping results 
in explosions that are spherically symmetric in CASTRO; 
low-order departures from spherical symmetry are sup- 
pressed. Our models therefore only capture higher-order 
asymmetries in the explosions. 

Perturbations that give rise to the Rayleigh- Taylor in- 
stabilities are seeded from the axisymmetric or Cartesian 
grid. We performed calculations in 2D at twice the res- 
olution of the simulations considered in this paper and 
compared the results of the high resolution simulations to 
the low resolution simulations. The width of the mixed 
region and velocities of the elements were essentially the 
same, indicating that our results are numerically con- 
verged and do not depend on the initial scale of pertur- 
bations imposed by the grid. 

3.3. Enlarging the Domain 

In order to minimize the amount of computational re- 
sources used, we implemented a strategy of enlarging the 
simulation domain only as necessary as the size of the re- 
gion of interest increased in time. The simulations were 
initialized on a 128™ grid with 2 levels of refinement, 
where n is the dimensionality of the simulation. This 
gives an effective grid resolution of 512 n at the finest 
level. A given simulation was advanced in time to the 
point where the shock was near the edge of the domain. 
The simulation was then stopped and the domain was 
doubled in each coordinate dimension. The enlarged do- 
main was covered by grids with cell spacing twice that 



4. RESULTS 

In each multidimensional simulation presented here, 
when the shock encountered the large region of increas- 
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ing pr 3 at the base of the hydrogen envelope of the 
presupernova star (or, in the case of model zl5, the 
helium-hydrogen envelope) it slowed, giving rise to a re- 
verse shock. The reverse shock then propagated inwards 
in mass towards the center of the star, and left a re- 
versed pressure gradient in its wake. This triggered the 
Ray leigh- Taylor instability first between the helium shell 
and the hydrogen envelope, and then between the helium 
and oxygen layers of the supernova. For model zl5, this 
instability first arose between the oxygen layer and the 
helium-hydrogen envelope. The instability grew, mixing 
together layers of increasing atomic mass as the reverse 
shock propagated inward through the star. Eventually 
the reverse shock passed by, and mixing stopped once the 
Raylcigh- Taylor fingers had dissipated their momenta. 

In each simulation, the domain was enlarged five times, 
bringing the final effective resolution of each simulation 
to 16, 384™, where n is the dimensionality of the sim- 
ulation. For model ul5, four levels of refinement were 
retained at the final stage of the simulation; for models 
sl5 and zl5, two levels of refinement were retained. 

Mixing had effectively ceased by 4.0xl0 5 seconds for 
model sl5, 3.8ccl0 4 seconds for model ul5, and 1.7xl0 5 
seconds for model zl5. Mixing had stopped by similar 
times in 2D and 3D simulations. The models were run to 
twice this time to make sure that no additional mixing 
would occur. 

Shown in Figure [5] are snapshots at identical times 
for the two- and three-dimensional simulations of models 
(from top) sl5, ul5, and zl5. The abundance of He, O, 
Si, and 56 Ni are shown in a logarithmic scale, starting 
clockwise from top left. These snapshots correspond to 
times when mixing had effectively stopped in the differ- 
ent models. On the left are snapshots of the 2D models; 
at right are snapshots showing slices in the YZ plane at 
the X axis. Upon first inspection, there appears to be 
very little difference between the two and three dimen- 
sional simulations of all three models. The width of the 
mixed region appears the same between the 2D and 3D 
models, although mixing in the 3D models appears to 
be more complete in this region than in the 2D mod- 
els. The shape of the instability differs slightly between 
the 2D and 3D simulations. The mushroom shape of 
the instability is more clearly defined in 2D than in 3D, 
especially for model ul5, where the instability had less 
time to grow and thus retained its original shape more 
clearly. In the 3D models, the Rayleigh- Taylor instability 
appears more elongated, in line with what was found in 
the single mode study of iKane et al.l (|2000l ). The great- 
est difference between the 2D and 3D simulations arises 
in model sl5, shown at the bottom of Figure [3J where 
there is actually less mixing of 56 Ni in 2D than in 3D. 

In all models the Rayleigh- Taylor fingers have grown 
to the point where they have begun to interact with 
one another. Models zl5 and sl5, in which the RT in- 
stability had the longest time to grow, show the great- 
est degree of interaction. This interaction transfers en- 
ergy and momentum in the transverse directions to the 
blast wave, leading to a transit ion to turbulence in 3D , 
and to a chaotic regime in 2D (jRemington et al.ll2006[ ). 
I Miles et all ()2005j ) state that this transition to turbulence 
begins to happen when the Rayleigh- Taylor fingers are 
about 5 to 6 times longer than their initial wavelength. 
This has clearly happened in all simulations. Models z!5 



and sl5 appear strikingl y similar to the fu lly turbulent 
simulations presented in iMiles et al.1 ()2005j ) . 

The differences between 2D and 3D can be compared 
more quantitatively by examining Figure [3] This figure 
displays the radial average of abundances for significant 
elements in both 2D and 3D simulations of the three 
models after all mixing has stopped. The times shown are 
the same times as in Figure[2j Figure[3]shows clearly that 
the width of the mixing region in the 2D and 3D models 
is nearly the same. 56 Ni is the exception, and is slightly 
more mixed in 2D than in 3D in the simulations of model 
sl5. Previous simulations of the Rayleigh- Taylor insta- 
bility in which these fingers do not interact show « 30% 
faster growth rate in 3D than in 2D, resulting in a wider 
mixed region in 3D. We also see a faster initial growth 
rate in 3D than in 2D, but once the Rayleigh- Taylor fin- 
gers begin to interact with one another, the growth rate 
in the 3D models decreases and the final width of the 
mixed region is the same between 2D and 3D. This is in 
line with simulations of the Rayleig h- Taylor ins t abilit y 
in a laboratory context presented in lMiles et al.l (|2005l ). 
in which the RT fingers interacted. 

The 2D and 3D simulations do appear somewhat differ- 
ent in Figure [3] in that the 3D lines are smoother than the 
2D lines. Because interactions between Rayleigh- Taylor 
fingers can transfer momentum and energy perpendicu- 
lar to the direction of the blast wave, mixing proceeds to 
a greater extent in 3D than in 2D. This reflects the fact 
that the 3D simulations have become fully turbulent (in 
the case of zl5 and sl5), while the 2D simulations are 
chaotic, as 3D simulations are able to transfer momen- 
tum and energy in two dimensions transverse to the blast 
wave, while 2D simulations have only one transverse di- 
mension at their disposal. There is likely a sampling 
effect at play here, as well. In 2D, there simply are not 
as many Rayleigh- Taylor fingers as there are in 3D, and 
thus there is effectively a smaller sample of the instabil- 
ities that grow. In 3D there are more Rayleigh- Taylor 
structures, so the space the fingers can occupy is better 
sampled, leading to a smoother radialized distribution. 

The velocities obtained when mixing has effectively 
ceased are virtually identical between the 2D and 3D 
simulations, as shown in Figure |U These simulations are 
shown at different times: 4.0a;10 5 seconds, 3.8a;10 4 sec- 
onds, and 1.7a;10 5 seconds for models sl5, ul5, and zl5, 
respectively, and so direct comparison of one model with 
another would be misleading, as the material will slow 
down as time goes on. It is clear from Figure H that 
the ultimate velocities after mixing has stopped and the 
Rayleigh- Taylor instability has frozen out are essentially 
the same i n the 2D and 3D simul ations. This is different 
from what lHammer et al.l (|2009l ) found in their simula- 
tions with larger degrees of initial asymmetry. They saw 
a high velocity tail for 56 Ni, O, and Si. We do not see a 
high velocity tail here. 

The 3D renderings shown in Figure [5] give a clearer 
picture of the shape of the Rayleigh- Taylor instabilities 
in our simulations than the slice through the 3D plane 
shown in Figure [5J We show renderings of the outermost 
surfaces of carbon, oxygen, silicon, and nickel. Model sl5 
is shown at the left, model ul5 in the center, and model 
zl5 at the right. These renderings were taken at different 
times, but after mixing had frozen out. Only in model 
z!5 have bits of the oxygen shell broken off and begun to 
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Fig. 2. — Snapshots of elemental abundance in the 2D (left) and 3D (right) versions of (from top to bottom) models sl5, ul5, and zl5. 
The 3D snapshots are slices at the X axis through the YZ plane. Because of artifacts in CASTRO, the extent of mixing is likely to be 
slightly exaggerated in slices at an axis. Nevertheless, the width of the mixed region in 2D and 3D appear very similar in all models. There 
is even slightly less mixing of 56 Ni in model sl5 in 3D than in 2D. The shape of the instability differs slightly between 2d and 3d in the 
expected manner-the 2D shapes arc slightly more "mushroomed" and the 3D shapes are more elongated, with the 3D models appearing 
to have transitioned to a fully turbulent regime, but otherwise the two simulations appear very similar. 
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Fig. 3. — Abundance by mass of individual elements as a function of mass coordinate. The width of the mixed region is essentially the 
same between the 2D and 3D cases, though 3D is smoother. The higher mass coordinate layers of the simulations have been omitted from 
the plot to better show the details of their inner portions. 
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Fig. 4. — Mass of individual elements as a function of velocity per 100 km s velocity bin. The amount of heavier elements such as O 
and 56 Ni at high velocities is essentially the same in 2D and 3D. Where the two differ, as they do most significantly for Ni and Si in models 
s!5 and z!5, the 2D simulations actually show material mixed out to higher velocities than the 3D simulations. 



move through the H-Hc shell that surrounds the heavy 
element core of this star. For no models have clumps 
of heavier elements penetrated shells of lighter elements, 
as happened i n the asymmetric explosion modeled by 
lHammer et al.l (|2009f ). Even the broken-off clumps of 
oxygen visible in the right hand panel of Figure [5] are 
encased in a layer of carbon. It is also apparent by ex- 
amining this figure that mixing has progressed further in 
the simulations where the progenitor star was a red giant, 
model sl5 and zl5, than for model ul5, whose progeni- 
tor died as a more compact blue star. The original shape 
of the instability can be seen in the carbon surface (top 
middle panel). In addition, the 56 Ni at the center of the 
model is less deformed in this model than in previous 
models. The most vigourous mixing has taken place in 
model zl5, which lacked a dense helium shell and died 
as a red giant. 

5. DISCUSSION 

Previous simulations of the Ray leigh- Taylor instability 
in core collapse supernovae have found significant dif- 
ferences between simil ar calculations per formed in two 
and three dimensions. iKane et al.l (|2000f ) examined the 
growth of Rayleigh- Taylor fingers arising from a single 
mode perturbation, and found that the instability grew 
30 — 35% faster in 3D than in 2D in a 1987A-type pro- 
genitor model. This faster growth rate was not sufficient 
to replicate obs ervations of 198 7 A, ho wever. In the more 
recent paper of lHammer et all (|2009() . the authors per- 
formed simulations of a 1987A-type model in two and 



three dime nsions, s tartin g with a more realistic initial 
model from lSch eck (2003). Their initial model was in the 
first stages of the explosion, and showed large-scale asym- 
metry in the heavy element core. They found, again, that 
single perturbations grew around 30% faster in 3D than 
in 2D. The authors proposed that artificial drag forces 
arising from 2D geometry slow down the growth of the 
instability in 2D simulations relative to 3D simulations, 
where these artificial drag forces a r e not present. This 
was also described in IKane et al.l (|2000D . who explain 
that the rising bubbles have less effective matter to push 
out of the way in 3D than i n 2D. 

The lHammer et ail (|2009t) study also found that in 2D, 
rising bubbles of heavy elements became entrained in the 
dense He shell left behind by the shock. In 3D, however, 
these bubbles were able to burst through the He shell to 
reach the H envelope, replicating one of the key observa- 
tions of 1987A. 

Studies of the Rayleigh- Taylor instability outside the 
context of supernova explosions have also found that a 
single Ray leigh- T aylor instability grows about 30% faster 
in 3D than in 2D ([Anuchina et al.ll2004D . Previous simu- 
lations o f the Rayleigh- Ta ylor instability in a laboratory 
context dMiles et al.ll2005l ) have found that the width of 
the mixed region is very similar between 3D and 2D, 
provided that the individual Rayleigh- Taylor fingers in- 
teract. Initially, the Rayleigh- Taylor fingers grow faster 
in 3D than in 2D. Once the fingers begin to interact, 
3D simulations transition to turbulence more completely 
than 2D simulations. 3D simulations can transfer mo- 
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Fig. 5. — 3D renderings of the outermost surface of carbon, oxygen, silicon, and nickel after mixing has frozen out. At left is model sl5; 
center is ul5, and right is zl5. Note that while the scale and time are consistent within models, the 3 models are shown at different scales 
and time: 5a;10 13 cm and A.OxlO 5 seconds for sl5, 5xl0 12 cm and 3.8xl0 4 seconds for model ul5, and 3xl0 13 cm and 1.7xl0 5 seconds for 
model zl5. For no model did heavier, inner layers penetrate the lighter, outer layers, though for model zl5 clumps of material have broken 
off and moved a small distance into the H-Hc envelope. 



10 



mentum in two dimensions perpendicular to the blast 
wave, while 2D simulations can transfer momentum in 
only one transverse direction. This results in more com- 
plete mixing within the mixed region in 3D than in 2D, 
reducing the Atwood number in 3D relative to 2D. 

For the simplest case of two incompressible fluids un- 
der constant acceleration with an initial perturbation of 
wavelength A and amplitude 77, such that 77 <C A, th e 
growth rate of the instability is given bv lRavleighl(fT88l : 



dt 



(7) 



where g is a constant acceleration and t is time. 
A, the Atwood number, is defined by 

A = ^ P ' 2 ~ 



(P2+Pl) 



(8) 



where p\ and P2 are the densities of the light and heavy 
fluids, respectively. 

As the densities of the two fluids become compera- 
ble, the Atwood number becomes smaller and the growth 
rate decreases. The growth rate of the bubbles in 3D is 
thus reduced, neatly canceling out the reduced drag these 
same bubbles experience in 3D as opposed to 2D. The 
width of the mixed region remains essentially the same 
between 2D and 3D for cases where the Ray leigh- Taylor 
fingers interact with one another. 

It could be supposed that the reason the final state of 
the supernova remnant appears the same in our 2D and 
3D models is that the perturbations arising from the grid 
were somehow larger in 2D than in 3D, and this effect 
exactly cancelled out the faster growth rates we expect 
to see in 3D. But this is not the case. Perturbations 
arising from the grid are larger in 3D than in 2D. These 
perturbations arise because we try to fabricate a spheri- 
cal object with boxes, as required by the RZ and Carte- 
sian geometries used in our simulations. The maximum 
amount one would expect to depart from a sphere (or cir- 
cle) in space is equivalent to the longest distance across a 
single cell, which in 2D is ^/2/S.x and in 3D is y/3Ax. To 
make certain, we performed the same set of calculations 
in 2D at twice the resolution. These simulations give the 
same results as our less refined simulations, indicating 
that our simulations are numerically converged for the 
problem presented here. 

As i n the si mulations oflMiles et al.l (|2005l ) , iKane et al.l 
(|200Q[) . and lHammer et all (|2009t ). we find that the 
Raylcigh- Taylor instability initially grows faster in 3D 
than in 2D. As the Raylcigh- Taylor fingers grow to the 
point where they begin to interact with one another the 
simulations begin to transition to a turbulent (3D) or 
chaotic (2D) regime. In 3D, momentum is transferred 
from the fingers in two dimensions, as opposed to only 
one dimension in 2D. The tranisiton to turbulence in 3D 
results in more complete mixing within the mixing region 
than the 2D simulation, reducing the Atwood number in 
3D relative to 2D. This has the effect of reducing the 
growth rate of the instability in 3D relative to 2D. The 
reduced growth rate in 3D relative to 2D at late times 
compensates for the fact that the bubbles do indeed ex- 
perience less drag in 3D than in 2D, and grow faster at 
early times. By the time mixing has ceased, the width 



of the mixed region is nearly identical in 2D and 3D in 
models where the Raylcigh- Taylor fingers have time to 
interact with one another. 

It can be seen from Eq. [7] that instabilitcs seeded by 
long wavelength perturbations grow more slowly than 
those seeded by short wavelength perturbations. Ob- 
viously, the interior of a supernova deviates from the 
ideal case: the acceleration is not constant, the fluids 
are compressible, and r\ soon becomes of the same or- 
der as A. In order for long wavelength perturbations 
to dominate before mixing ceases there must be signif- 
icant initial power at long wavelengths. Because there 
is less time for the instability to grow in blue stars, we 
expect to see more traces of large perturbations at fairly 
low wavenumb ers in these sta r s than in red stars. The 
simulations of Ha mmer et all (|2009f) have a significant 
amount of power at long wavelengths, so these pertur- 
bations grow fast enough to dominate before mixing is 
truncated. 

The initial perturbations in our simulations are a great 
deal smaller than those examined in previous studies. 
The spectrum of these perturbations is also flatter and 
made up of higher modes. The perturbations are more 
random, though, it should be noted, not entirely so, and 
are largest around the 30 and 60° marks on the grid. This 
can be seen during the initial stages of the development 
of the instability, though not at later times, for model 
si 5. The enhanced perturbations at these points is likely 
behind the placement of the clumps visible in the 3D 
rendering of model zl5 in Figure [5] 

Mixing had more time to develop in our two red gi- 
ant stars, models zl5 and sl5, than it generally does 
in the blue giant 198 7A-ty pe progenitor models studied 
in IKane et all (|2000l) and lHammer et all (|2009f ). Mix- 
ing froze out after 4xl0 5 and 1.7x10 s seconds for the 
red models sl5 and zl5, respectively, o ver an order of 
magni tude later than in the simulation of lHammer et all 
(|2009l) . in which mixing froze out after 9xl0 3 seconds. 
In our blue ul5 model, mixing froze out after 4xl0 4 sec- 
onds, which is still n early 5 times longer than in the 
lHammer et all ([2009D model. The individual Rayleigh- 
Taylor instabilities in all simulations grew for sufficient 
times that they interacted, though more interaction took 
place in the red models. This was not the case in th e 
single-mode perturbation explored in Ka ne et all (l2000f) . 
or th e few large perturbation modes in (Ha mmer et all 
120091 ). Smaller wavelengths grow faster than longer 
wavelengths, so longer wavelength perturbations are less 
likely to become nonlinear. In these models, as well 
as the single instability investigated in lAnuchina et all 
(|2004[) , the Ray leigh- Taylor fingers did not interact with 
one another, ensuring that the relative drag on the fingers 
was the only significant difference between 2D and 3D. 
This accounts for the around 30% faster growth rates 
discovered in these simulations. The simplicity of the 
initial perturbation spectrum can delay the transition to 
turbulence, and he nce make mixing appear more efficient 
in 3D than in 2D (| Miles et al.ll2005H 

For long enough wavelength perturbations, spherical 
geometry will ensure that the perturbations diverge, so 
they never interact and thus never enter a turbulent 
regeime. This will quickly lead to significant depar- 
tures from axisymmetry. There may exist a transition 
mode perturbation which would be the longest wave- 
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length mode that could lead to turbulent interaction. 
The instabilities will diverge without interaction if the 
wavelength of an instability is greater than the height at 
which it would become nonlinear. This tranistion mode 
ought to be dependent on presupernova structure: red 
supcrgiants have longer time for the Raylcigh- Taylor in- 
stability to grow, and thus potentially a lower mode than 
blue giants. 

The Rayleigh- Taylor instability may "forget" its initial 
conditions at late times. This can happen though an in- 
verse cascade, where shorter wavelength bubbles merge 
to form longer wavelength structures. Under idealized 
conditions, the flow may become self similar. This tran- 
sition to self-similarity likely happens during the early 
stages of Rayleigh- Taylor growth in our simulations, ac- 
counting for the similar end appearance of our simula- 
tions despite increasing the resolution (and decreasing 
the initial perturbations). This inverse cascade can only 
lead to the instabilities becoming so big. A necessary 
condion for self-similarity is that 

A <C /i < L (9) 

where h is the height of the instability and L is the 
size of the region in which the instability is growing. 
The merger of short wavelength instabilities to larger 
wavelength instabilites should be truncated long before 
A becomes comperablc to L, which would produce large- 
scale asymmetries. Large-scale departures from symme- 
try must arise from the supernova explosion mechanism 
itself, as they cannot arise through an inverse cascade of 
the Rayleigh- Taylor instability. 

Recent observations suggest that asymmetry in 
core-collapse type sup ernova explosions is common 
(|Wang k Wheeled 12008( 1. For SNe type II, this asym- 
metry is characterized by a dominant polarization an- 
gle, which is most likely due to a directed, jet-like 
flow. Polarization data show that most SNe also exhibit 
composition-dependent, large scale departures from ax- 
isymmctry. The degree of departure from a xisy mmetry 
may b e less in SNell-P: IChugai et all (|2005l > and lChugail 
(120011 are able to match observations of one SNell- 
P, SN2004dj, by using a bipolar distribution of 56 Ni. 
This fits with the emerging picture of collapse supernova 
mechanisms, which are lik ely to be inherentl y mulAtidi- 
mensional and aspher ical (iJanka et al.ll2007f) The SASI 
explosion mechanism dBurrows et al.l 120061 : IScheck et al.l 
120061 iMarek k Jankall2009| l seems to be one of the better 
candidates for reviving the stalled shock and successfully 
exploding a massive star, and it is inherently asymmet- 
ric. How relevant, then, are spherical explosion models? 
This is the first time the post-explosion hydrodynamics 
of a wide range of supernova progenitors has been stud- 
ied in 3D. It is reasonable to model Pop III and Pop II 
core collapse SNe in 3D first as spherical explosions as 
no observations of these primordial supcrnovae have yet 
been made. This study will contextualize future studies 
of supcrnovae progenitors with asymmetric explosions. 

For type II SNe in which departures from asymmetry 
are small and dominated by an axisymmetric component, 
2D simulations may do as well as 3D simulations in repro- 
ducing the axisymmetric features of the supernova rem- 
nant. In these stars, Rayleigh- Taylor fingers of all but 
the lowest mode order are expected to grow far enough 



that they begin to interact with one another. Simulations 
in 2D of such stars can be used to predict nucleosynthe- 
sis and fallback for suc h explosions with a fa ir degree of 
accuracy, as was done in lJoggerst et al.l (|2010T ) . 2D calcu- 
lations may also be used to more accurately predict light 
curves of SNe for which departures from axisymmetry do 
not dominate, which would use far less computational re- 
sources than similar calculations performed in 3D. 

6. CONCLUSIONS 

We have performed two and three dimensional simu- 
lations of supernova models with three different mctal- 
licities: model zl5, with Z = Z@, model ul5, with 
Z= 10~ 4 Zq, and model sl5, with Z=Zq. The initial 
perturbations are seeded by the grid. Calculations per- 
formed at twice the resolution in 2D indicate that our 
solutions are numerically converged. We find no signifi- 
cant difference in the width of the mixed region between 
our 3D and 2D simulations of the same supernova. Pre- 
vious studies of Raylcigh- Taylor mixing in 1987A-type 
supernova models found that mixing was more efficient 
and individual Raylcigh- Taylo r instabilitie s grew about 
30% faster in 3D than in 2D. iKane et al.l (|2000t l stud- 
ied single-mode density perturbations with amplitudes of 
10% of mode I = 10, and found that perturbations grew 
3 0-35% faster in 3D th an in 2D. The more recent study 
of lHammer et al.l (|2009l ) examined the post-explosion hy- 
drodynamics of a 1 987A-type mo del with a fully 3D ex- 
plosion model from lScheckl (|2007f ) with large scale asym- 
metries at early times. They found that the instability 
grew about 30% faster in 3D than in 2D. 

While we also find a faster initial growth rate in 3D 
than in 2D in our simulations, the growth rate in 3D de- 
clines once the simulations transition to turbulence, and 
the final width of the mixed region is the same in 2D 
and 3D. Because the Raylcigh- Taylor instability in our 
model is seeded by perturbations arising from the grid, 
our perturbations have a lower amplitude, and a flatter 
spectrum shifted tow ards higher modes tha n either the 
IKane et al.l (|2000t l or lHammer et~all (|2009tl studies. In 
our progenitor models the Raylcigh- Taylor instability has 
more time to develop and grow than in previous simula- 
tions. With a growth time on the order of 10 5 seconds, 
instabilities in models sl5 and zl5 grow for more than 
10 tim es longer than the simulations of lHammer et all 
(|2009H . in which mixing froze out at around 9x10 sec- 
onds. In model ul5, which died as a blue supcrgiant, 
mixing ceases after 4xl0 4 seconds, which is a shorter time 
than the red giant models in this survey b ut is still nearly 
5 time s longer than the the simulations of lHammer et all 
<[2009h . 

In our progenitor models and with perturbations 
seeded by the grid the RT instability has enough time 
to grow that the fingers of the instability begin to inter- 
act with one another. This interaction transfers momen- 
tum and energy transverse to the shock, allowing the 3D 
simulations to mix more fully than 2D simulations. The 
more efficient transfer of momentum and energy in 3D 
renders the 3D simulations fully turbulent in the mixed 
region. This higher degree of mixing in 3D than in 2D 
reduces the Atwood number in the 3D simulations rel- 
ative to the Atwood number in the 2D simulations. A 
lower Atwood number leads to a reduction in the bub- 
ble growth rate in 3D, which compensates for the lower 
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effective drag in 3D. The result is that the width of the 
mixed region is virtually identical between simulations 
performed in 2D and those performed in 3D. 3D simula- 
tions have transitioned to full turbulence and are more 
completely mixed than 2D simulations. Thi s was found 
in a la boratory context in the simulations of iMilcs et al.1 
(|2005f ). For simulations in which interactions between 
the Rayleigh- Taylor fingers do not occur, either b ecause 
only one finger is modeled (|Anuchina et al.ll2004D or the 
Rayleigh- Taylor fing ers grow for a shor t enough time that 
they never interact (jKane et al.ll2000T) . the greater effec- 
tive drag experienced in 2D versus 3D will continue to 
dominate the simulation, ensuring that the width of the 
mixed region will be larger in 3D than in 2D. For simula- 
tions in which the Rayleigh- Taylor instability has enough 
time to grow that the fingers begin to interact with one 
another, the width of the mixed region will be the same 
in 2D and 3D simulations. 

Instabilities arising from long wavelength perturba- 
tions take longer to grow than those arising from short 
wavelength perturbations. Large amplitude perturba- 
tions of long wavelength are necessary for such instabili- 
ties to grow substantially. These large-scale fingers may 
diverge before they can transition to turbulence, and thus 
may give rise to departures from axisymmetry. The large 
amplitude, long wavelength pertu rbations present in the 
more realistic explosion model of lHammer et all (|2009fl 
ensure that such fingers don't interact, which results in 
the significant differences between their 2D and 3D sim- 
ulations. Rayleigh- Taylor theory predicts that an inverse 
cascade, in which shorter wavelength instabilities merge 
to form longer wavelength instabilitcs, should be trun- 
cated long before the wavelcgth A becomes comparable 
to the size of the region in which the instability is grow- 
ing. Large scale departures from symmetry do not arise 
through an inverse cascade in our 3D simulations; they 



must have their origins in the explosion itself. 

It is likely that many, if not most, core-collapse super- 
novae explosions are highly asymmetric. Observations 
not just of 1987A but of other supernova, as well the 
"kicks" observed in pulsars, indicate that asymmetry is 
a near-universal phenomenon amongst the explosions of 
massive stars. From the theoretical side, it seems in- 
creasingly likely that viable supernova explosions, like 
the SASI mechanism, are inherently asymmetric. In su- 
pernovae that die as red giants the Rayleigh- Taylor in- 
stability will have enough time to develop that individual 
fingers will interact with one another. If the asymmetries 
in such explosions are dominated by axisymmetric com- 
ponents, 2D models might well be good predictors of the 
width of the mixed region and the axisymmetric features 
of the remnant. Future work should investigate asym- 
metric explosions in a wide variety of progenitor models. 
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